packs.inst <- c("readxl","foreign","dplyr","tidyr","ggplot2","stargazer","haven",
                "dummies","Hmisc","lmtest","sandwich", "doBy", "multiwayvcov",
                "miceadds", "car", "purrr", "knitr", "zoo", "readstata13", 
                "tidyverse", "psych","wesanderson", "lubridate","reporttools",
                "data.table", "devtools","rmarkdown","estimatr", "ivpack", 
                "Jmisc", "lfe", "plm", "tinytex", "xts", "psych", 
                "PerformanceAnalytics","roll", "rollRegres", "glmnet", "hdm", 
                "broom", "RCurl", "learnr", "maps", "fGarch", "remotes", 
                "RPostgreSQL", "DBI", "RPostgreSQL", "remotes", "RPostgres", 
                "Rmisc", "ggthemes", "splitstackshape", "gginference", 
                "MASS","fBasics","xtable")
packs.load <- c("fGarch")
# lapply(packs.inst, install.packages, character.only = FALSE) 
lapply(packs.inst, require, character.only = TRUE)

Part 1: deal with the data

import data

# real inflation rate from monthly year-on-year CPI
CPI <- read.csv("~/Desktop/DTFF_final/data/CPI.csv")
# expected inflation rate from 1-year bond yield to maturity
bond <- read.csv("~/Desktop/DTFF_final/data/bond.csv")
# commodity futures' price data
commodity <- read.csv("~/Desktop/DTFF_final/data/commodity.csv")[-1,]
# stock and gold price
stock_gold <- read.csv("~/Desktop/DTFF_final/data/stock_gold.csv")[-1,]
# real estate price
realestate <- read.csv("~/Desktop/DTFF_final/data/real_estate.csv")

time-series data

dateymd <- as.Date(ymd(commodity[,1]))
CPI_ts <- xts(x = CPI[,-1], order.by = dateymd)
bond_ts <- xts(x = bond[,-1], order.by = dateymd)
commodity_ts <- xts(x = commodity[,-1], order.by = dateymd)
stock_gold_ts <- xts(x = stock_gold[,-1], order.by = dateymd)
realestate_ts <- xts(x = realestate[,-1], order.by = dateymd[18:156])

merger the returns

# calculate the monthly year-on-year logarithm yield
colnames(CPI_ts) <- "monthly year-on-year CPI"
CPI_log_ts <- log(1+CPI_ts[,1])
colnames(CPI_log_ts) <- "real_inflation"
# 1-year bond t-1 the yield of maturity is used as the 
# expected inflation rate at time t
bond_lag_ts <- stats::lag(bond_ts,1)
colnames(bond_lag_ts) <- "expected_inflation"
# monthly year-on-year logarithm yield of commodities
commodity_log_ts <- log(commodity_ts) - log(stats::lag(commodity_ts,12))
# monthly year-on-year logarithm yield of stock and gold
stock_gold_log_ts <- log(stock_gold_ts) - log(stats::lag(stock_gold_ts,12))
# monthly year-on-year logarithm yield of real estate
realestate_log_ts <- log(realestate_ts) - log(stats::lag(realestate_ts,12))

CPI_log_ts <- CPI_log_ts[13:156,]
bond_lag_ts <- bond_lag_ts[13:156,]
commodity_log_ts <- commodity_log_ts[13:156,]
stock_gold_log_ts <- stock_gold_log_ts[13:156,]
unexpected_inflation <- CPI_log_ts - bond_lag_ts 
colnames(unexpected_inflation) <- "unexpected_inflation"
realestate_log_ts <- realestate_log_ts[13:139,]
DATA3 <- merge.xts(CPI_log_ts, bond_lag_ts, unexpected_inflation, 
                   commodity_log_ts, stock_gold_log_ts, realestate_log_ts)
save(DATA3,file = "DATA.RData")

Part 2: Examination of the Inflation Hedging Ability

regression of commodity futures

## commodity futures
library(stargazer)
cm_lm1 =list()
for (i in 4:7) {
  cm_lm1[[i-3]]  <- lm(DATA3[,i] ~ DATA3$expected_inflation + DATA3$unexpected_inflation)}
# cm_reg1 <- stargazer(cm_lm1, dep.var.labels = "commodity futures", align = TRUE, df = FALSE)
stargazer(cm_lm1, dep.var.labels = "commodity futures", align = TRUE, df = FALSE, 
          type = "text")
## 
## =========================================================
##                              Dependent variable:         
##                      ------------------------------------
##                               commodity futures          
##                         (1)       (2)      (3)     (4)   
## ---------------------------------------------------------
## expected_inflation   -5.983*** -8.515*** 4.295*   -1.201 
##                       (2.063)   (2.253)  (2.480) (2.207) 
##                                                          
## unexpected_inflation   0.280    1.704*    0.928  -2.303**
##                       (0.876)   (0.957)  (1.053) (0.937) 
##                                                          
## Constant             0.205***  0.251***  -0.078   0.012  
##                       (0.057)   (0.062)  (0.068) (0.061) 
##                                                          
## ---------------------------------------------------------
## Observations            144       144      144     144   
## R2                     0.071     0.157    0.021   0.042  
## Adjusted R2            0.058     0.145    0.007   0.029  
## Residual Std. Error    0.138     0.151    0.166   0.148  
## F Statistic          5.366***  13.119***  1.530  3.103** 
## =========================================================
## Note:                         *p<0.1; **p<0.05; ***p<0.01
cm_lm2 =list()
for (i in 8:11) {
  cm_lm2[[i-7]]  <- lm(DATA3[,i] ~ DATA3$expected_inflation + DATA3$unexpected_inflation)}
# cm_reg2 <- stargazer(cm_lm2, dep.var.labels = "commodity futures", align = TRUE, df = FALSE)
stargazer(cm_lm2, dep.var.labels = "commodity futures", align = TRUE, df = FALSE,
          type = "text")
## 
## ==========================================================
##                               Dependent variable:         
##                      -------------------------------------
##                                commodity futures          
##                        (1)      (2)       (3)       (4)   
## ----------------------------------------------------------
## expected_inflation   -2.744  -9.813*** -10.418*** -5.093**
##                      (2.186)  (2.878)   (2.420)   (2.575) 
##                                                           
## unexpected_inflation  0.616   2.774**   2.056**   -2.176**
##                      (0.928)  (1.222)   (1.028)   (1.042) 
##                                                           
## Constant              0.094  0.302***   0.311***  0.152** 
##                      (0.060)  (0.079)   (0.067)   (0.072) 
##                                                           
## ----------------------------------------------------------
## Observations           144      144       144       140   
## R2                    0.021    0.157     0.193     0.043  
## Adjusted R2           0.007    0.145     0.182     0.029  
## Residual Std. Error   0.146    0.193     0.162     0.164  
## F Statistic           1.540  13.171*** 16.912***   3.049* 
## ==========================================================
## Note:                          *p<0.1; **p<0.05; ***p<0.01
cm_lm3 =list()
for (i in 12:15) {
  cm_lm3[[i-11]]  <- lm(DATA3[,i] ~ DATA3$expected_inflation + DATA3$unexpected_inflation)}
# cm_reg3 <- stargazer(cm_lm3, dep.var.labels = "commodity futures", align = TRUE, df = FALSE)
stargazer(cm_lm3, dep.var.labels = "commodity futures", align = TRUE, df = FALSE,
          type = "text")
## 
## ==========================================================
##                               Dependent variable:         
##                      -------------------------------------
##                                commodity futures          
##                        (1)      (2)       (3)       (4)   
## ----------------------------------------------------------
## expected_inflation   -7.728** -4.968**  -4.505   -7.502***
##                      (3.181)  (2.188)   (2.936)   (1.733) 
##                                                           
## unexpected_inflation  -1.291  -1.685*  -5.235*** 5.103*** 
##                      (1.351)  (0.929)   (1.247)   (0.736) 
##                                                           
## Constant             0.251*** 0.157***  0.145*   0.269*** 
##                      (0.087)  (0.060)   (0.081)   (0.048) 
##                                                           
## ----------------------------------------------------------
## Observations           144      144       144       144   
## R2                    0.040    0.042     0.111     0.425  
## Adjusted R2           0.027    0.029     0.098     0.417  
## Residual Std. Error   0.213    0.146     0.196     0.116  
## F Statistic           2.952*  3.113**  8.806***  52.096***
## ==========================================================
## Note:                          *p<0.1; **p<0.05; ***p<0.01
cm_lm4 =list()
for (i in 16:19) {
  cm_lm4[[i-15]]  <- lm(DATA3[,i] ~ DATA3$expected_inflation + DATA3$unexpected_inflation)}
# cm_reg4 <- stargazer(cm_lm4, dep.var.labels = "commodity futures", align = TRUE, df = FALSE)
stargazer(cm_lm4, dep.var.labels = "commodity futures", align = TRUE, df = FALSE,
          type = "text")
## 
## =======================================================
##                             Dependent variable:        
##                      ----------------------------------
##                              commodity futures         
##                         (1)       (2)     (3)     (4)  
## -------------------------------------------------------
## expected_inflation   -14.338*** -4.748  -2.444   0.377 
##                       (3.966)   (4.215) (3.704) (3.346)
##                                                        
## unexpected_inflation  3.739**   -1.961  -1.508  -0.351 
##                       (1.685)   (1.790) (1.535) (1.387)
##                                                        
## Constant              0.392***   0.095   0.080   0.015 
##                       (0.109)   (0.116) (0.102) (0.093)
##                                                        
## -------------------------------------------------------
## Observations            144       144     142     140  
## R2                     0.166     0.013   0.008   0.001 
## Adjusted R2            0.154    -0.001  -0.007  -0.014 
## Residual Std. Error    0.265     0.282   0.241   0.218 
## F Statistic          14.007***   0.897   0.533   0.057 
## =======================================================
## Note:                       *p<0.1; **p<0.05; ***p<0.01

regression of spot gold

## spot gold
sg <- lm(DATA3$Gold_price ~ DATA3$expected_inflation + 
           DATA3$unexpected_inflation)
# sg_reg <- stargazer(sg, dep.var.labels = "Spot Gold", align = TRUE, df = FALSE)
stargazer(sg, dep.var.labels = "Spot Gold", align = TRUE, df = FALSE,
          type = "text")
## 
## ================================================
##                          Dependent variable:    
##                      ---------------------------
##                               Spot Gold         
## ------------------------------------------------
## expected_inflation            -7.454***         
##                                (1.714)          
##                                                 
## unexpected_inflation          5.128***          
##                                (0.728)          
##                                                 
## Constant                      0.267***          
##                                (0.047)          
##                                                 
## ------------------------------------------------
## Observations                     144            
## R2                              0.431           
## Adjusted R2                     0.423           
## Residual Std. Error             0.115           
## F Statistic                   53.374***         
## ================================================
## Note:                *p<0.1; **p<0.05; ***p<0.01

regression of industry stocks

## industry stocks
stock_lm1 =list()
for (i in 21:25) {
  stock_lm1[[i-20]]  <- lm(DATA3[,i] ~ DATA3$expected_inflation + DATA3$unexpected_inflation)}
# stock_reg1 <- stargazer(stock_lm1, dep.var.labels = "industry stocks", align = TRUE, df = FALSE)
stargazer(stock_lm1, dep.var.labels = "industry stocks", align = TRUE, 
          df = FALSE, type = "text")
## 
## ===================================================================
##                                   Dependent variable:              
##                      ----------------------------------------------
##                                     industry stocks                
##                        (1)      (2)       (3)       (4)      (5)   
## -------------------------------------------------------------------
## expected_inflation   -2.188  -9.063**   -7.416    -3.864    -4.955 
##                      (3.515)  (4.328)   (4.488)   (3.589)  (3.551) 
##                                                                    
## unexpected_inflation -2.023  -6.328*** -6.728*** -5.326***  -1.565 
##                      (1.493)  (1.838)   (1.907)   (1.525)  (1.509) 
##                                                                    
## Constant             -0.017   0.250**   0.206*    0.169*   0.300***
##                      (0.097)  (0.119)   (0.123)   (0.099)  (0.098) 
##                                                                    
## -------------------------------------------------------------------
## Observations           144      144       144       144      144   
## R2                    0.013    0.082     0.082     0.080    0.016  
## Adjusted R2          -0.001    0.069     0.069     0.067    0.002  
## Residual Std. Error   0.235    0.290     0.300     0.240    0.238  
## F Statistic           0.925  6.296***  6.287***  6.134***   1.126  
## ===================================================================
## Note:                                   *p<0.1; **p<0.05; ***p<0.01
stock_lm2 =list()
for (i in 26:30) {
  stock_lm2[[i-25]]  <- lm(DATA3[,i] ~ DATA3$expected_inflation + DATA3$unexpected_inflation)}
# stock_reg2 <- stargazer(stock_lm2, dep.var.labels = "industry stocks", align = TRUE, df = FALSE)
stargazer(stock_lm2, dep.var.labels = "industry stocks", align = TRUE, 
          df = FALSE, type = "text")
## 
## ===================================================================
##                                   Dependent variable:              
##                      ----------------------------------------------
##                                     industry stocks                
##                         (1)        (2)      (3)     (4)      (5)   
## -------------------------------------------------------------------
## expected_inflation   -11.363***  -2.650   -6.993*  4.802   -4.069  
##                       (3.302)    (3.117)  (4.168) (4.398)  (3.401) 
##                                                                    
## unexpected_inflation -4.228***  -5.384*** -2.201  -1.412  -7.689***
##                       (1.403)    (1.324)  (1.770) (1.868)  (1.445) 
##                                                                    
## Constant              0.416***    0.090   0.253** -0.135    0.085  
##                       (0.091)    (0.086)  (0.115) (0.121)  (0.094) 
##                                                                    
## -------------------------------------------------------------------
## Observations            144        144      144     144      144   
## R2                     0.098      0.108    0.023   0.019    0.171  
## Adjusted R2            0.085      0.095    0.009   0.005    0.159  
## Residual Std. Error    0.221      0.208    0.279   0.294    0.228  
## F Statistic           7.641***  8.541***   1.625   1.390  14.542***
## ===================================================================
## Note:                                   *p<0.1; **p<0.05; ***p<0.01

regression of real estate

## real estate
estate_lm = list()
for (i in 31:33) {
  estate_lm[[i-30]]  <- lm(DATA3[,i] ~ DATA3$expected_inflation + DATA3$unexpected_inflation)}
# estate_reg <- stargazer(estate_lm, dep.var.labels = "real estate", align = TRUE, df = FALSE)
stargazer(estate_lm, dep.var.labels = "real estate", align = TRUE, 
          df = FALSE, type = "text")
## 
## ==================================================
##                           Dependent variable:     
##                      -----------------------------
##                               real estate         
##                         (1)       (2)       (3)   
## --------------------------------------------------
## expected_inflation     0.720    2.110**  2.589*** 
##                       (1.294)   (0.850)   (0.835) 
##                                                   
## unexpected_inflation  -0.172     0.449    0.704*  
##                       (0.568)   (0.373)   (0.367) 
##                                                   
## Constant               0.046    -0.010    -0.030  
##                       (0.036)   (0.024)   (0.023) 
##                                                   
## --------------------------------------------------
## Observations            127       127       127   
## R2                     0.004     0.049     0.079  
## Adjusted R2           -0.012     0.034     0.064  
## Residual Std. Error    0.077     0.051     0.050  
## F Statistic            0.277    3.189**  5.338*** 
## ==================================================
## Note:                  *p<0.1; **p<0.05; ***p<0.01

result of inflation hedging effect

assetnames <- colnames(DATA3[,-(1:3)])
Expected   = c(assetnames[c(1,2,3,6,7,8,9,10,12,13,17,19,23,25,29,30)],NA)
Unexpected = assetnames[c(2,4,6,7,8,10,11,12,13,17,19,20,21,23,24,27,30)]
HedgingAbility <- cbind(Expected, Unexpected)
as.data.frame(HedgingAbility)
xtable(as.data.frame(HedgingAbility))

Part 3: Inflation Hedgig Portfolio by Mean Variance Method

exclude assets which have no inflation hedging effect

asset_returns <- DATA3[,-c(1:3)]
asset_returns <- asset_returns[,-c(5,14,15,16,18,22,26,28)]
save(asset_returns, file="asset_returns.RData")
dim(asset_returns)
## [1] 144  22

the set of attainable portfolios

# Calculate the mean vector and covariance matrix
mean_ret <- colMeans(asset_returns, na.rm = T)
cvar_ret <- cov(na.omit(asset_returns))

# Calculate individual weights 
for (i in colnames(asset_returns)[-(dim(asset_returns)[2])]) {
  weight = runif(10000, min=-1.5, max=1.5)
  names = paste0(i,"_weight")
  if (i == "soybeans.No..1"){
    weight_final = weight
    names_final = names}
  else {
    weight_final = cbind(weight_final, weight)
    names_final = cbind(names_final, names)}
  }

# Get the dataframe and matrix on the weights
weight_df <- as.data.frame(weight_final)
colnames(weight_df) <- names_final

weight_df$sum <- rowSums(weight_df)
weight_df$third.tier <- 1 - weight_df$sum 
weight_df$sum <- NULL

matrix_weights <- as.matrix(weight_df)

# report the weights of the matrix
str(data.frame(matrix_weights))
## 'data.frame':    10000 obs. of  22 variables:
##  $ soybeans.No..1_weight               : num  1.399 0.508 -0.668 -1.257 0.317 ...
##  $ soybeans.No..2_weight               : num  1.012 0.585 0.257 -0.195 -1.368 ...
##  $ yellow.corn_weight                  : num  0.528 -0.21 1.394 1.407 0.992 ...
##  $ LLDPE_weight                        : num  1.3392 1.1248 0.9379 -1.0823 0.0362 ...
##  $ palm.oil_weight                     : num  -0.5 -0.58 -0.42 -1.435 -0.986 ...
##  $ soybean.oil_weight                  : num  1.4 1.215 -0.373 -1.197 -0.973 ...
##  $ PVC_weight                          : num  -0.738 0.374 -0.424 0.705 0.738 ...
##  $ cathode.copper_weight               : num  -1.394 1.033 -0.504 -0.414 1.001 ...
##  $ aluminum_weight                     : num  0.964 0.776 1.28 0.963 -1.379 ...
##  $ zinc_weight                         : num  1.305 0.932 -0.329 0.352 -0.732 ...
##  $ gold_weight                         : num  -0.0944 0.4378 -0.56 -0.4133 1.1655 ...
##  $ natural.rubber_weight               : num  1.1774 -0.4654 0.9295 0.0138 -0.605 ...
##  $ Gold_price_weight                   : num  1.414 1.352 -0.219 -0.197 -1.03 ...
##  $ Material_price_weight               : num  0.979 0.33 1.446 -1.077 -1.17 ...
##  $ Industry_price_weight               : num  -0.151 1.019 0.82 -1.382 -0.794 ...
##  $ Optional_consumption_price_weight   : num  -0.149 -1.238 0.321 -0.165 0.356 ...
##  $ Medicine_health_care_price_weight   : num  1.282 0.825 -0.33 -1.239 -0.132 ...
##  $ Finance_and_real_estate_price_weight: num  1.36 -1.468 0.491 1.487 -1.4 ...
##  $ Information_technology_price_weight : num  -0.8762 0.5551 0.0913 0.5529 -0.3066 ...
##  $ Utility_price_weight                : num  -0.385 -1.289 -1.166 1.144 0.242 ...
##  $ second.tier_weight                  : num  -1.477 -0.343 -0.647 -1.014 -1.143 ...
##  $ third.tier                          : num  -7.4 -4.47 -1.33 5.44 8.17 ...
head(data.frame(matrix_weights))
# Calculate the feasible expected returns and standard deviations
feasible_pf_mu = matrix_weights%*%mean_ret
feasible_pf_sd = apply(matrix_weights, 1, 
                       function(x) sqrt(t(x) %*% cvar_ret %*% x))

# Construct the feasible dataframe
# consisting of 10000 differently weighted risk and return combinations
feasible_pf <- as.data.frame(cbind(feasible_pf_mu, feasible_pf_sd))  
colnames(feasible_pf) <- c("Portfolio_Return", "Portfolio_Risk")

# report the feasible dataframe
str(feasible_pf)
## 'data.frame':    10000 obs. of  2 variables:
##  $ Portfolio_Return: num  -0.083843 0.038258 -0.076509 -0.000588 0.160481 ...
##  $ Portfolio_Risk  : num  1.123 0.88 0.778 0.952 1.145 ...
head(feasible_pf)
# Now, let's visualise the relationship
feasible_pf %>% 
  ggplot(aes(x= Portfolio_Risk, y = Portfolio_Return)) + 
  geom_point(color = "grey") + 
  geom_point(data = subset(feasible_pf, Portfolio_Risk <= 0.12 & 
                             Portfolio_Return >= 0), color = "darkorchid3", 
             shape = 1, aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.12 & 
                             Portfolio_Risk <= 0.14 & Portfolio_Return >= 0.07), 
             color = "darkorchid3", shape = 1,aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.14 & 
                             Portfolio_Risk <= 0.16 & Portfolio_Return >= 0.09), 
             color = "darkorchid3", shape = 1,aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.16 & 
                             Portfolio_Risk <= 0.18 & Portfolio_Return >= 0.1), 
             color = "darkorchid3", shape = 1,aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.18 & Portfolio_Risk 
                           <= 0.20 & Portfolio_Return >= 0.11), 
             color = "darkorchid3",shape = 1, aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.2 & 
                             Portfolio_Risk <= 0.22 & Portfolio_Return >= 0.11), 
             color = "darkorchid3", shape = 1, aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  ylab(expression(mu[p])) + xlab(expression(sigma[p])) + 
  ggtitle("Porfolio Frontier for 10'000 different weights") +
  labs(color='Factor Portfolios') +
  theme(plot.title= element_text(size=14, color="grey26",
  hjust=0.3,lineheight=2.4, margin=margin(15,0,15,0)), 
  panel.background = element_rect(fill="#f7f7f7"),
  panel.grid.major.y = element_line(size = 0.5,linetype ="solid",color="grey"),
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  plot.background = element_rect(fill="#f7f7f7", color = "#f7f7f7"), 
  axis.title.y = element_text(color="grey26", size=12, margin=margin(0,10,0,10)),
  axis.title.x = element_text(color="grey26", size=12, margin=margin(10,0,10,0)),
  axis.line = element_line(color = "grey"))

## the minmum - variance portfolio

# Define the matrix A. It consists of: 

## the covariance matrix multiplied by two
## a column right to the covariance matrix, consisting of 1's
## a row right below the covariance matrix and the additional column, 
## consisting of 1's and one zero (the zero is in the right-bottom 
## of the resulting matrix)
mat_A <- rbind(cbind(2*cvar_ret, rep(1, dim(cvar_ret)[1])), 
               c(rep(1, dim(cvar_ret)[1]), 0))
# Define the vector b as vector of zeros with dimension of the covariance 
# matrix and one 1 at the bottom
vec_b <- c(rep(0, dim(cvar_ret)[1]), 1)
# Calculate the inverse and perform matrix multiplication to get the vector z
z <- solve(mat_A)%*%vec_b
# Derive the first N elements of the vector to retrieve the actual values 
x_MV <- z[1:dim(cvar_ret)[1]]
# Check that the sum adds up to 1
sum(x_MV)
## [1] 1
# Calculate the expected return: 
mu_MV <- x_MV %*% mean_ret
sd_MV <- sqrt(t(x_MV) %*% cvar_ret %*% x_MV)
# Create the appropriate dataframe
MV_PF <- as.data.frame(cbind(mu_MV, sd_MV, t(x_MV)))
colnames(MV_PF) <- c("Mu_MV", "Sd_MV",names_final, "third.tier_weight")
as.data.frame(t(MV_PF))
# Now, let's visualize the relationship
feasible_pf %>% 
  ggplot(aes(x= Portfolio_Risk, y = Portfolio_Return)) + 
  geom_point(color = "grey") +
  # This is just to color in the "optimal PFs"
  geom_point(data = subset(feasible_pf, Portfolio_Risk <= 0.12 & 
                             Portfolio_Return >= 0.02), color = "darkorchid3", 
             shape = 1, aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.12 & 
                             Portfolio_Risk <= 0.14 & Portfolio_Return >= 0.07), 
             color = "darkorchid3", shape = 1,aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.14 & 
                             Portfolio_Risk <= 0.16 & Portfolio_Return >= 0.09), 
             color = "darkorchid3", shape = 1,aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.16 & 
                             Portfolio_Risk <= 0.18 & Portfolio_Return >= 0.1), 
             color = "darkorchid3", shape = 1,aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.18 & 
                             Portfolio_Risk <= 0.20 & Portfolio_Return >= 0.11), 
             color = "darkorchid3",shape = 1, aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  geom_point(data = subset(feasible_pf, Portfolio_Risk > 0.2 & 
                             Portfolio_Risk <= 0.22 & Portfolio_Return >= 0.11), 
             color = "darkorchid3", shape = 1, aes(x= Portfolio_Risk, y = Portfolio_Return)) +
  # Calculate and plot the Minimum Variance PF
  geom_point(color = "goldenrod",  aes(x= MV_PF$Sd_MV, y = MV_PF$Mu_MV), 
             size = 3) +
  annotate('text',x = -0.01 ,y = 0, label = "Minimum Variance PF", 
           size = 3.5, color = "goldenrod") +
  ylab(expression(mu[p])) + xlab(expression(sigma[p])) + 
  ggtitle("Minimum Variance portfolio") +
  labs(color='Factor Portfolios') +
  xlim(-0.10, 0.75) + 
  theme(plot.title= element_text(size=14, color="grey26",
  hjust=0.43,lineheight=2.4, margin=margin(15,0,15,0)), 
  panel.background = element_rect(fill="#f7f7f7"),
  panel.grid.major.y = element_line(size = 0.5, linetype = "solid",color="grey"),
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  plot.background = element_rect(fill="#f7f7f7", color = "#f7f7f7"), 
  axis.title.y = element_text(color="grey26", size=12, margin=margin(0,10,0,10)),
  axis.title.x = element_text(color="grey26", size=12, margin=margin(10,0,10,0)),
  axis.line = element_line(color = "grey"))
## Warning: Removed 4848 rows containing missing values (geom_point).

## portfolio efficient frontier

# First we calculate the first Efficient Portfolio
# Define the EW return
mu_spec_x <- mu_MV
mu_spec_y <- 0.1

# We first define again the Matrix A 
mat_A_EF <- rbind(cbind(2*cvar_ret, mean_ret,  rep(1,dim(cvar_ret)[1])), 
                  cbind(t(mean_ret), 0, 0), 
                  cbind(t(rep(1,dim(cvar_ret)[1])), 0, 0))

# Then, we define the vector b
vec_b_EF_x <- c(rep(0, dim(cvar_ret)[1]), mu_spec_x, 1)
vec_b_EF_y <- c(rep(0, dim(cvar_ret)[1]), mu_spec_y, 1)

# Now, we can solve for the respective weights
z_EF_x <- solve(mat_A_EF)%*%vec_b_EF_x
z_EF_y <- solve(mat_A_EF)%*%vec_b_EF_y

# Then, we take the first N elements to get the respective weights
x_EF <- z_EF_x[1:dim(cvar_ret)[1]]
y_EF <- z_EF_y[1:dim(cvar_ret)[1]]

# Now, let's calculate the risk
sd_EF_x <- sqrt(t(x_EF) %*% cvar_ret %*% x_EF)
sd_EF_y <- sqrt(t(y_EF) %*% cvar_ret %*% y_EF)
sd_EF_xy <- t(x_EF) %*% cvar_ret %*% y_EF
# Lastly, compute the weights and results
a = seq(from=0, to=1, by=0.05) 
# Create the expected return as well as the variance and standard 
# deviation of each portfolios
mu_z = a * mu_spec_x + (1-a) * mu_spec_y
## Warning in a * mu_spec_x: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
sd_z = sqrt(a^2*sd_EF_x^2 + (1-a)^2*sd_EF_y^2 + 2*a*(1-a)*sd_EF_xy)
## Warning in a^2 * sd_EF_x^2: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning in (1 - a)^2 * sd_EF_y^2: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning in 2 * a * (1 - a) * sd_EF_xy: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
# Create a dataframe consisting of different weights
z = matrix(0, length(a), dim(asset_returns)[2]) 

for (i in 1:length(a)){
  z[i, ] = a[i]*x_EF + (1-a[i])*y_EF
  }

# Create a dataframe consisting of only the efficient linear transformation 
# portfolios
z_df <- as.data.frame(cbind(mu_z, sd_z))
colnames(z_df) <- c("Efficient_PF_Return", "Efficient_PF_Risk")
# Now, let's visualise the relationship
z_df %>% 
  ggplot(aes(x= Efficient_PF_Risk, y = Efficient_PF_Return), color = "goldenrod") + 
  geom_point() + 
  geom_path() + 
  ylab(expression(mu[p])) + xlab(expression(sigma[p])) + 
  ggtitle("The Portfolio Efficient Frontier") +
  labs(color='Factor Portfolios') +
  theme(plot.title= element_text(size=14, color="grey26",
  hjust=0.3,lineheight=2.4, margin=margin(15,0,15,0)), 
  panel.background = element_rect(fill="#f7f7f7"),
  panel.grid.major.y = element_line(size = 0.5, linetype="solid", color="grey"),
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  plot.background = element_rect(fill="#f7f7f7", color = "#f7f7f7"), 
  axis.title.y = element_text(color="grey26", size=12, margin=margin(0,10,0,10)),
  axis.title.x = element_text(color="grey26", size=12, margin=margin(10,0,10,0)),
  axis.line = element_line(color = "grey")) 

# Calculate the TP
## HINT: Choose a similar risk free rate as in the exercise session. 
### import the risk-free data set
rf <- DATA3[,2]
### calculate the TP
mu_f = mean(rf)
## First the numerator and denominator 
numerator_T_N <- inv(cvar_ret) %*% (mean_ret - mu_f*rep(1, length(mean_ret)))
denominator_T_N <- t(rep(1, length(mean_ret))) %*% 
  inv(cvar_ret) %*%(mean_ret - mu_f*rep(1, length(mean_ret)))
## calculate the weights
weights_T_N <- numerator_T_N[,1] / denominator_T_N
## Warning in numerator_T_N[, 1]/denominator_T_N: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
weights_T_N
##                soybeans.No..1                soybeans.No..2 
##                   0.092002203                  -0.563613294 
##                   yellow.corn                         LLDPE 
##                  -0.018099874                  -0.468700215 
##                      palm.oil                   soybean.oil 
##                   0.005661339                   0.445684816 
##                           PVC                cathode.copper 
##                  -0.376695617                   0.438619433 
##                      aluminum                          zinc 
##                   0.947920112                  -0.335028039 
##                          gold                natural.rubber 
##                   0.877784704                  -0.594614747 
##                    Gold_price                Material_price 
##                  -0.564797082                  -0.447550157 
##                Industry_price    Optional_consumption_price 
##                   0.666089223                   0.464717306 
##    Medicine_health_care_price Finance_and_real_estate_price 
##                   0.675056516                   0.203704284 
##  Information_technology_price                 Utility_price 
##                  -0.813196211                  -0.625157650 
##                   second.tier                    third.tier 
##                   4.440973314                  -3.450760364
return_T_N <- weights_T_N %*% mean_ret
sd_T_N <- sqrt(t(weights_T_N) %*% cvar_ret %*% weights_T_N)
# Create the tangent portfolio
### define the risk of the risk-free asset and the covariance 
### of it with risky assets
sigma_f = 0
sigma_Af = 0

## Define the sequence
x_tan = seq(from=-0.8, to=2.2, by=0.1)
x_f = 1 - x_tan

## Calculate the metrics
### Calculate the risk and return metrics
mu_T_N <- weights_T_N %*% mean_ret
sd_T_N <- sqrt(t(weights_T_N)%*%cvar_ret%*%weights_T_N)

### Create another dataframe
mu_sd_T_N_df <- as.data.frame(cbind(mu_T_N, sd_T_N))
colnames(mu_sd_T_N_df) <- c("mu_T_N", "sd_T_N")

mu_tanf = x_tan*mu_T_N + x_f*mu_f
## Warning in x_tan * mu_T_N: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
var_tanf = x_tan^2*sd_T_N^2 + x_f^2*sigma_f^2 + 2*x_tan*x_f*sigma_Af
## Warning in x_tan^2 * sd_T_N^2: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
sd_tanf = sqrt(var_tanf)

# Only get the "positive returns" from the r_f on
mu_tanf_real <- mu_tanf[1:9]
sd_tanf_real <- sd_tanf[1:9]
# Create a df
cml_N <- as.data.frame(cbind(mu_tanf_real, sd_tanf_real))
colnames(cml_N) <- c("MU_CML", "SD_CML")
mu_sd_T_N_df
port <- as.data.frame(cbind(Minimum_Variance = c(t(x_MV), mu_MV, sd_MV),
                            Tangency = c(weights_T_N, return_T_N, sd_T_N)))
rownames(port) <-c(colnames(asset_returns), "Return", "sd")
xtable(port, caption = "Weights of two portfolios", digits = 5)